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Abstract 

Spurred by an experimental controversy in the literature, we investigate the end-monomer dynamics 
of semiflexible polymers through Brownian hydrodynamic simulations and dynamic mean-field theory. 
Precise experimental observations over the last few years of end-monomer dynamics in the diffusion 
of double-stranded DNA have given conflicting results: one study indicated an unexpected Rouse-like 
scaling of the mean squared displacement (MSD) {r^{t)) ~ t^^^ at intermediate times, corresponding to 
fluctuations at length scales larger than the persistence length but smaUer than the coil size; another study 
claimed the more conventional Zimm scaling {r^{t)) ~ t^^^ in the same time range. Using hydrodynamic 
simulations, analytical and scaling theories, we find a novel intermediate dynamical regime where the 
effective local exponent of the end-monomer MSD, a{t) = dlog(r^(t))/dlogt, drops below the Zimm 
value of 2/3 for sufficiently long chains. The deviation from the Zimm prediction increases with chain 
length, though it does not reach the Rouse limit of 1/2. The qualitative features of this intermediate 
regime, found in simulations and in an improved mean-field theory for semiflexible polymers, in particular 
the variation of a{t) with chain and persistence lengths, can be reproduced through a heuristic scaling 
argument. Anomalously low values of the effective exponent a axe explained by hydrodynamic effects 
related to the slow crossover from dynamics on length scales smaller than the persistence length to 
dynamics on larger length scales. 

1 Introduction 

Recent experimental advances using fluorescence correlation spectroscopy^"^ have given unprecedented in- 
formation about the dynamical behavior of large single polymer molecules in solution, in particular the 
small-scale kinetics of individual monomers inaccessible to traditional techniques like dynamic light scatter- 
ing. One of the first studies along this direction yielded an unexpected result. Shusterman et al? observed 
the random motion of a single labeled monomer at the end of a long double-stranded DNA molecule, and 
found evidence of an "intermediate Rouse regime" : the mean squared displacement (MSD) followed a scaling 
{r^{t)) oc t^/'^ for a wide time range corresponding to polymer motion at length scales smaller than the coil 
size Rg but larger than the persistence length Ip. This agrees with the free-draining Rouse model for a poly- 
mer which neglects hydrodynamic interactions mediated by flow fields arising from the monomers moving 
through the solvent. Such a result contradicts the conventional wisdom for flexible polymers, which states 
that these hydrodynamic interactions play a crucial role in polymer dynamics in dilute solutions and give rise 
to non-draining behavior that is qualitatively described by the Zimm theory, which predicts (r^{t)) cx t^/^.^'"^ 
Though double-stranded DNA is a semiflexible polymer (having a persistence length w 50 — 100 nm much 
larger than the width « 2 nm), the expectation for kinetics at scales larger than Ip is that it behaves like a 
non-draining flexible polymer. Thus the apparent absence of hydrodynamic effects is quite surprising, and 
the intermediate Rouse regime does not fit into established theories of the dynamics of flexible polymers in 
dilute solutions, though recently there has been an attempt to explain its existence through a theory exhibit- 
ing time-dependent hydrodynamic screening.^ On the other hand, a new experimental study by Petrov et. 
al.^ on the same system did not seem to show the Rouse regime, and its results were interpreted to be gen- 
erally consistent with the dynamics predicted by the Zimm theory. Arguably the dynamics of a semiflexible 
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polymer such as DNA may be expected to differ from tfiat of flexible polymers. However, with the exception 
of the Harnau, Winkler, Reineker (HWR) model, ^ other established theories of the dynamics of semi-flexible 
polymers^°' treat only the range of displacements smaller than the persistence length {r'^{t)) < Ip. 

To help resolve the controversy over the dynamics of semiflexible polymers on intermediate length scales 
Ip < {r^{t)) < L^, we study the end- monomer behavior of semiflexible chains in dilute solutions using 
two approaches: dynamic mean-field theory (MFT) that includes hydrodynamics with the pre-averaging 
approximation, and Brownian hydrodynamics simulations without the pre-averaging approximation. The 
end- monomer MSD, diffusion constants, and longest relaxation times from the two approaches agree closely 
with each other. While the hydrodynamic pre-averaging MFT method is similar to that of HWR in Ref. 9, 
we have improved the approximation by taking into account the full hydrodynamic interaction matrix in the 
Langevin equation, and not just the diagonal contribution. This leads to much better agreement between the 
MFT and the simulation data: compared with the earlier version, the improved MFT is 10-65% closer to the 
mean-square displacement (r^(i)) of the end monomer obtained by simulations for time scales shorter than 
the longest relaxation time of the chain, and reduces the discrepancy in the effective local exponent a{t) = 
d\og{r^ (t)) / dlogt in this time range, which is underestimated by as much as 10% using the earlier method. 
Thus we can confidently extend the MFT to larger chain lengths L that are inaccessible to simulation. 
For these chains we find an intermediate dynamical regime where the continuously varying effective local 
exponent of the end-monomer MSD, a{t), drops below 2/3, and its difference from this Zimm value increases 
with L. The existence of this regime and the qualitative trends of a(t) with changing L and Ip are verified 
independently through a heuristic scaling argument. However even at the largest chain lengths examined, 
comparable to or longer than the experimentally studied chains of Refs. 2 and 5, the effective exponent a(t) 
does not reach the Rouse limit of 1/2. Comparison with the experimental MSD data of Ref. 2 reveals two 
interesting results: the MFT accurately describes the long-time diffusion behavior, related to the large-scale 
dynamics of the chain; however at shorter times it underestimates the extent of the MSD. As we show in 
this paper, the same sub-Zimm scaling of the MSD is also contained in the HWR theory that was used 
to successfully fit the data of Ref. 5. So the question is not whether an intermediate sub-Zimm scaling 
regime exists, but rather how large that regime is and how small are the intermediate exponents. The 
remaining discrepancy between theory and experiment discussed in this paper highlights the importance of 
additional dynamical degrees of freedom absent in the worm-like chain model used as the starting point for 
the theoretical description of DNA or some shortcomings in the current analysis of FCS data. 

The paper is organized as follows: in Sec. [5] we give a heuristic scaling argument that captures the basic 
properties of the intermediate dynamical regime; in Sec. [3] we describe the details of the Brownian dynamics 
simulations; in Sec. 3] we give an overview of the mean-field model for semiflexible polymers and the pre- 
averaging approximation used to determine its behavior in solution; in Sec. [5] we compare the simulation, 
MFT, and heuristic results, together with the experimental data. The dynamical regimes exhibited in these 
results are examined through asymptotic scaling analysis in Sec. [6l and placed in the context of earlier 
theories. Finally Sec. [7| summarizes the main points of the paper. Additional material, extending the mean- 
field model of Sec. |3]to extensible worm- like chains, is provided in Appendix A. Mathematical details of an 
analytical approximation used in Sec. [6] are given in Appendix B. 

2 Heuristic Scaling Argument 

Certain qualitative features of the intermediate dynamical regime for semifiexible polymers can be derived 
from a scaling argument similar to the one that is typically used to understand subdiffusive motion in the 
Rouse or Zimm models. In our heuristic scaling, we assume that at time scales t the diffusion of the 
polymer is characterized by the coherent motion of a section with contour length £{t). The MSD of a 
monomer during this time t has two limiting subdiffusive behaviors depending on the magnitude of £(t). 
For £{t) <C Ip, where the local stiffness of the polymer plays the key role, the MSD {r^{t)) is dominated by 
transverse chain fiuctuations and thus {r'^{t)) ^ £{t)^ /Ip}'^^^^ For i{t) S> Ip, the MSD can be estimated as 
{r^it)) ^ r1^{£{t)), where r1^{t} is the mean squared end-to-end distance of a polymer of contour length L 
For a semifiexible chain with persistence length Ip this is given by: 

rl{l)^2lpl-2ll{l-e-'/'^). (1) 
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Figure 1: Results of the heuristic scahng argument described in Sec. ^ (a) The effective exponent a = 
d\og{r^ (t)) / dlogt versus t/r for polymers of diameter 2a, contour length L = 12000a, and various persistence 
lengths Ip = 50a — 4000a. The dot along each a{t) curve marks a{Tp). The time scales r and Tp are defined in 
the text, (b) The minimum value of a over all t versus Ip/ L for several L = 200a— 12000a. (c) The minimum 
value of a over all t versus L/lp for several Ip = 50a — 1800a. (d) The circles show the minimum value of a 
in each of the curves in panel (b), plotted as a function of l/ln(L/a). Superimposed is a straight-line fit to 
the data points, min/^ mint ct « 0.48 + 1.04/ ln(i/a). 
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Figure 2: The contours show the minimum value of the effective exponent a{t) = dlog{r'^(t)) /dlogt over 
ah t calculated using the heuristic scaling argument in Sec. [5] for chains of diameter 2a, total length L, and 
persistence length Ip. The interval between contours is 0.02. The dotted line indicates the minimum value 
of mint ct over Ip/a for a given L/a. 



Thus the total intramolecular contribution to the MSD, 
smoothly interpolates between these two limits: 



j^(i)), can be written in a heuristic form which 



{rUm)y 



-1/01 



(2) 



Using the fact that r1^{t) ~ i"^ for £ <^ Ip, one can easily check that Eq. ^ has the appropriate asymptotic 
limits for small and large £{t), with the crossover occurring at t « Tp where the relaxation time Tp of a 



persistent segment is defined by r^ci^irp)) 



;2 



We added numerical constants Ci, C[ ^ 0(1) and a 



crossover exponent which will be determined from the comparison with numerical data. 

In order to close the scaling argument, we write the diffusion relation between spatial and temporal scales 



6D{e{t))t, 



(3) 



in terms of an effective diffusion constant D{£) for a polymer section of length £. D{£) can be estimated 
from a heuristic formula^® that compares well with both simulation and experimental DNA results: 



D{£)=C,iDt;,i£)+Dt^,,{£)r/'>'^ 



(4) 



where 



^cyl(^) 



^coil(^) 



2^ofc_Br 

£ 



la 



0.312 



+0.565 
2.22^0^5^0 



0.1 



(5) 



lp£ 



Here 2a is the diameter of the chain, = l/(67r?7a) is the Stokes mobility of a single sphere of radius a, 
and rj is the viscosity of water. Dcyi{£) is the diffusion constant of a stiff cylinder,^'' and Z^coiK-^) gives the 
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diffusion constant of an ideal chain in the flexible limit where £ ^ Ip. We assume £ <^ Ip/O''^, so self-avoidance 
effects can be ignored. For parameter values a — 1 nm and Ip = 50 nm, typical of DNA, these effects 
are only important at scales greater than ^ 10^ /ini, far larger than the chain lengths investigated in the 
experimental studies discussed above. Eq. Q is an interpolation between the limiting cases given by Dcyi{£) 
when £ <^ Ip and I3coii(-^) when £ ^ Ip. In Ref. 16 the crossover exponent was determined to be 02 = 3 and 
with the definition of Dcy\{£) and i?coii(^) the constant was fixed at C2 = 1. In the present context we use 
4>2 and C2 ^ 0(1) as parameters that will be adjusted to fit the numerical data. 

The asymptotic Zimm scaling is easily obtained from the expressions written so far: In the flexible regime, 
for £ > Ip, we have Dcoii ~ (^^p)~^^^ and from Eqs. ^ and ^ we find (r-^j^^j^) ~ £lp. Using (r^j^j.^^) ~ Dt from 
Eqs. ^ we thus obtain £ ~ 1?!'^ jlp and {r\^^^^ ^ t^l'^ , the well-known Zimm scaling for flexible polymers. "'^^ 
In the stiff polymer regime, for £ <^ Ip, we have Dcyi ^ log{£/a) and from Eq. ^ we find (r^^tra) ~ ^^/^p- 
Again using (ri^t^^) ~ Dt we obtain this time £ ~ {tiplogty/* and thus (rf^tTii) ~ t^^'^lp^^'^ log^^"^ t. The 
MSD scaling of semiflexible polymers has pronounced logarithmic corrections in the stiff polymer regime due 
to hydrodynamic effects. 

To get an expression for the total MSD, {r^{t)), that includes the long-time regime, one must consider 
also the crossover which occurs near time t = t, where £{t) ~ L, the total contour length of the chain. 
For t ^ T the effective diffusion constant is D{L), and {r'^{t)) ~ 6D{L)t, describing the trivial diffusion of 
the whole polymer coil, t corresponds approximately to the longest relaxation time of the polymer. This 
crossover is captured by yet another crossover expression, 

(r'it)) = [{rLrM)^' + {6D{L)t}'^^] '^^^ , (6) 

which gives the correct asymptotic scaling behavior for {r'^{t)) in all time regimes. In the results below, the 
exponents (j)2, '/'a, and the three constants Ci, C[, C2 in Eqs. ([2]), (|4]), and ([6]) are chosen so that the 
heuristic scaling argument approximately agrees, for long chain lengths, with the numerical results of the 
MFT approach described in the next section. The best-fit values are: 4>i = 0.84, 02 = 3.15, ips = 3.62, 
C\ = 3.65, C[ = 0.66, C2 = 1.31. 

The full time dependence covering also the non-asymptotic behavior is obtained by equating the expres- 
sions for in Eqs. ^ and ^ and implicitly solving for £{t); we thus calculate {rf^tiai^)) ^ function 
of t. Plugging the result for (^^^^^.^^(t)) into Eq. ^ gives the total MSD {r^{t)). The time evolution of {r'^{t)) 
can be expressed through the effective exponent a{t) = d\og{r^{t)) /dlogt. Fig. [T]^a) shows a{t) versus t/r 
for chains of total length L = 12000a, with various persistence lengths in the range Ip = 25a — 4000a. The dot 
along each a(t) curve marks a(Tp). There is clearly an intermediate time regime, within the range Tp < t < t, 
where a{t) dips below the Zimm value of 2/3. The minimum value of a{t) over all t depends both on L and 
Ip, as shown in Figs.^b) and (c), which plot mint o; versus Ip/L for several chain lengths L — 100a— 12000a, 
and mint ct versus L/lp for several Ip in the range Ip = 50a — 1800a. The overall variation of min^ a as a 
function of L /a and Ip/ a is depicted in the contour diagram of Fig. (5] The deviation from Zimm behavior 
becomes more prominent with increasing L: the time range where a < 2/3 increases, and the values of min^ a 
decrease. As seen in Fig. [TJc), for fixed Ip/a the decrease in mint ct with L eventually saturates for L ^ Ip. 
The mint ct curves in Fig. [TJb) all reach a minimum in the range Ip/L ~ 0.01 — 0.04. The position of the 
minimum decreases with L approximately with the logarithmic dependence Ip/L w —0.013 -I- 0.26/ log(L/a). 
The exponent values at these minima, min; mint a, also have a nearly linear dependence on l/log(L/a), 
as can be seen in Fig. [U^d), where minj miuf a goes from 0.677 ai L — 200a to 0.591 at L = 12000a. The 
best-fit line is min;^ mint ce « 0.48-1-1.04/ log(L/a). (If data from L much larger than the experimental range 
is also included, the L = 00 extrapolation of minj^ mint a shifts from 0.48, approaching 1/2.) The growing 
deviation from Zimm behavior with L is possibly related to the observation in Ref. 2 that the intermediate 
Rouse regime becomes more noticeable at longer coil sizes, occupying a larger range of times. However, in 
contrast to Ref. 2, the exponent a never reaches the true Rouse value of 1/2 even at the longest realistic 
chain lengths. For Ip — 50a, corresponding to the DNA persistence length, mint a at L = 12000a is 0.602. 

The origin of this intermediate regime where a{t) < 2/3 can be linked to the crossover behaviors of 
(''Ltra) ^^'^ Assume £ is in the range Ip < £ < L and is sufficiently large that {rf^t,.^) ~ £lp, but 

small enough that D{£) has not reached the asymptotic limit £'coii(^) ^ i£lp)~^^'^- Since the total MSD 
(r^) ^ (^fntra) ^^^^ regime, one can use (r^^tra) ~ -^^ ^'^^ (^fntra) ^^p relate the effective exponent 
a — d\og{r'^) / dlogt to D{£), giving a = (1 — {£/ D)dD / d£)~^ . As D{() is in the crossover region between 
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Dcyi(^) and Z)coii(^), it must decrease with I slower than (. ^\og{i/a), but faster than i ^1"^. These two 
Hmits mean that a is bounded by 1/(2 - log"^(^/a)) > 1/2 from below, and 2/3 fr om above, corresponding 
precisely to the intermediate dynamical regime. 

The existence of this regime will be confirmed through the Brownian dynamics and MFT calculations 
described in the next two sections. In Sec. [S] we will see that the qualitative trends illustrated in Fig. [TJ^a)- 
(d) agree very well with the results from the more sophisticated MFT approach and thus allow for a simple 
explanation of sub-Zimm scaling behavior in terms of hydrodynamic effects on the diffusion behavior of a 
semiflexible polymer in the crossover between two limiting regimes. 



3 Brownian Dynamics Simulation 

For the numerical Brownian dynamics simulations^^' we model the polymer as a connected chain of M 
spheres, each having radius a and position Vi(t\i — 1, . . . , M. The sphere positions evolve in time according 
to the Langevin equation. 



dt 



M 

i=i 



dU{vi, . . .,vm) 



at): 



(7) 



appropriate for the low Reynolds number regime. Here "jji ij is the Rotne-Prager tensor^^ describing hydro- 
dynamic interactions between the monomers, 




(8) 



where Vij = — r^ , and 1 is the 3x3 identity matrix. The stochastic velocity ^j(t) in Eq. ([7]) is Gaussian, 
with correlations given by the fluctuation-dissipation equation: 

{Ut) ® ijit')) = ^kBT^jr.jdit - 1') . (9) 

The final component of the model is the elastic potential C/(r i , . . . , r j\/ ) in Eq. ([7]) , depending on the positions 
of the spheres. This potential consists of two parts. 



U ^ C/wLC + C/lj : 



(10) 



with 



M-l 



UwLC = ^ X! (''^+1:* ^ 2a)^ + ^ X! ~ cos 6*,;) 



[/lj =cj^e(2a-ry) 



2a 

i=2 

2ay^ /2a^^ 
r-ij J \ r,. 



(11) 



Here 9i is the angle between r^+i^i and ri_i_i. The C/wlc term describes the stretching and bending forces 
associated with the extensible worm-like chain model, with stretching modulus 7 and bending modulus e. 
The latter is related to the persistence length Ip of the polymer through e — IpkBT. For all the simulations 
the stretching modulus is set at 7 = 200kBT/a, which is large enough that the total contour length of the 
polymer stays approximately constant. The C/lj term is a truncated Lennard- Jones interaction with strength 
(jj = iksT. 

To implement Eq. ([7]) numerically, we discretize it with time step r, and use non-dimensionalized variables, 
measuring lengths in units of a, times in units of /{ksT^o), and energies in units of ksT. For a given 
contour length L = 2aM and persistence length Ip, the results described below are based on averages taken 
from 15 — 50 independent runs, each with time step t — 3 x 10^** / [ksT ^q) and lasting for ~ 10* — 10^ 
steps. The first 10^ steps of a run are not used for data collection, and afterwards output data are collected 
every 10^ — 10*^ steps. 
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4 Mean-field Model of Polymer Dynamics 



The derivation of the mean-field model for semiflexible polymers is described below. Readers not interested 
in the technical details may skip this section, the main result of which is Eq. ([57)1 for the end- monomer MSD 
in terms of several parameters: the diffusion constant D, relaxation times t„ and coefficients A„. All of 
these parameters can be determined for a given L and Ip by obtaining the normal modes and numerically 
evaluating the hydrodynamic interaction matrix H as outlined in Eqs. pop - (|34p . 

The analytical model of the polymer is a continuous space curve r(s) of total length L, with contour 
coordinate s in the range —L/2 < s < L/2. The simplest expression for the elastic energy U of the chain, 
incorporating the effects of rigidity, is that of Kratky and Porod,^^ 



2 J \ ds ^ 

where e is the bending modulus introduced above, e — IpksT, and the tangent vector u = dr/ds is subject 
to the constraint u^(s) = 1 at each s. As in Sec. Owe assume L ^ Ip/a'^, so we ignore self- avoidance effects. 
The associated free energy is = — /3^^ logZ, with f3~^ = ksT and the partition function Z given by the 
functional integral, 

Z= fvul[S{u^s)-l)e-''^ . (13) 

S 

The delta function enforcing the constraint can be equivalently written using an additional functional integral 
over a complex auxiliary field A(s), 



'too 



iOQ 



(14) 



too 



where we introduce the functional .F[A], and ignore any constants arising from the normalization of the 
integral. Since calculations with this partition function are generally intractable due to the tangent vec- 
tor constraint, we employ the mean-field theory (MFT) approach developed by Ha and Thirumalai,^'^' 
evaluating the functional integral over A(s) using a stationary-phase approximation: 

Z= / I?Ae-^^[^l«e"'^^[^-l. (15) 



Here Aci(s) is the path satisfying the stationary-phase condition 5T/5\{s) = 0, and we have neglected 
higher-order correction terms. The resulting MFT free energy Fmf = -^[^^ci] takes the form:^"^ 

Fmf = log J Pue^^^"^ , (16) 

where 



+ iyo{u\L/2) + u\~L/2)) , 



and the constants e, j/, and vq are related by: 

Y = ^o = ^fci3T. (18) 

Comparing the results of this mean-field analytical model to those of the simulation described in the last 
section, we note that the simulation potential energy in Eqs. (|10 |) - pT|) contains an additional extensional 
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term with large parameter 7 ^ ksT / a. Applying the mean-field approach to an extensible worm-like chain 
leads to a value of the effective stretching moduli v and vq slightly modified from that of Eq. p8)) . with 
corrections of the order of e/F, where T = 4a^7. The details are described in Appendix A. Assuming e 
and r are fixed in the continuum limit, and F 3> e, one can ignore the finite extensibility of the chain in 
constructing the mean-field theory. 

The MPT elastic energy of Eq. ([T7]) can be derived in several alternative ways: it was first proposed 
by Lagowski, Noolandi, and Nickel^^ as a modification of the Harris-Hearst model^^ that corrected chain 
inhomogeneities due to end fluctuations; it was later independently derived from the maximum entropy 
principle by Winkler, Reineker, and Harnau.^^ The main consequence of the approximation is that the local 
constraint u^(s) = 1 is relaxed and replaced by the condition (u^(s)) =1. If the relationship between the 
bending modulus e and Ip is redefined as e = {3/2)lpkBT in Umf, the tangent vector correlation function 
has the same form as in the Kratky-Porod chain. 



(u(s) • u(s')) = exp(-|s' - s\/lp) , 



(19) 



Related quantities like the mean squared end-to-end distance and radius of gyration are also correctly 
reproduced by the MPT approximation with this redefinition of e, and thus we will use it for the remainder 
of the paper. This applies only to the MPT elastic energy of Eq. p7|) : in the simulation Uwlc of Eq- (UHl) e 
retains its original definition. 

In deriving the diffusion behavior of the polymer in solution, we follow an approach similar to that 
of HWR,^ who first studied the dynamical characteristics of the MPT model given by Eq. P?|) using a 
hydrodynamic pre-averaging approximation along the lines of the Zimm model. ^'^ To describe the time 
evolution of the chain r(s,t) in the presence of hydrodynamic interactions, we start with the Langevin 
equation: 



di 



r{s,t) 



L/2 



ds' ^ (s, s'; r(s, t) — r(s', t)) 



-L/2 



(20) 



MF 



dr{s',t) 



Here the ^(s, t) is the stochastic contribution, and ^(s, s'; x) is the continuum version of the Rotne-Prager 
tensor in Eq. ([8]),^ 



^(s,s';x) = 2afioS{s — s') 1 
-|-e(a;- 2a) 



( ' 


^ X (g) X 
1 1 


\87r?7a; 




? 


x(Ki 


X 









(21) 



with the O function excluding unphysical configurations. 

The pre-averaging approximation consists of replacing jj, (s, s'; r(s, i) — r(s', t)) in Eq. (|20p . which involves 
a complicated dependence on the specific chain configuration at time t, with an average over all equilibrium 
configurations, /Xavg(s,s'), that depends only on the contour coordinates s and s' . This tensor /x avg is 
defined as: 



^avg(s, s') = j d^x ^?(s, s'; x)G'(s, s'; x) , 



(22) 



where G(s, s'; x) is the equilibrium probability of finding two points at s and s' along the polymer contour 
whose spatial positions differ by the vector x. Por the MPT model of Eq. (fT7|) . this probability is:^'^ 



where a{l) = 2lpl~2lp{l~exp{—l/lp)), the mean squared end-to-end distance of a chain of length I. Plugging 
Eq. ((231) into Eq. (HH wc find: 



M avg(Sj S ) — 



2a^QS{s — s') 



&{\s-s'\-2a) 
77 y^67rV(|s^^s^ 



' exp 



6a2 



(24) 



a{\s-s'\) 

=Mavg(s - s') 1 ■ 

For the same reason as in Eq. ([2T|) . we have added a Q function to the final result. 
The pre-averaged version of the Langevin equation is thus 

9 , ^ 



(25) 



V2 \ <5r(s',t). 

We assume the ^(s,t) are Gaussian random vectors, whose components f*^*^(s,i) have correlations given by 
the fluctuation-dissipation theorem: 

(^«(s,t)e(^)(s',0> =2fcBT5„-^(t-t')A'avg(s-s'). (26) 

Using J7mf from Eq. (fT7|) . the force term in Eq. ([25]) can be written as 

(5?7mf 



(5r(s',t) 

with free-end boundary conditions at s = ±i /2 of the form 



54 92 
^5^r(s',t)-)-2^— r(.',i) 



53 9 
e^r(±L/2, <) - 2^.— r(±L/2, i) , 

Te|^r(±i/2, i) - 2z.o|-r(±L/2, i) = . 
as^ as 



(27) 



(28) 



To rewrite the Langevin equation in matrix form, we assume ^(s,t) satisfies similar boundary conditions 
to r(s,i), and expand both r(s,i) and ^(s,t) in normal modes ipn{s), with amplitudes Pn{t) and qn(i) 
respectively: 



(29) 



We choose the normal modes V'n(s) to be eigenfunctions of the differential operator in Eq. ((27|) . satisfying 
for eigenvalues A„. These V'n(s) take the form:^ 



(30) 



■00 (s) = 

tpnis) = 



1 

L ' 

An I sm g^s 
L \ "cosQ!„L/2 
sinh/3„s 



cosh/3„L/2 



. odd, 



(31) 



cos a„s 



+f3n- 



sin a„L/2 
cosh /3„s 



sinh/3„L/2y ' 



n even. 
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0:^ - al = 2vlt , Ao = , A„ = eal + 2vai . (32) 



with 

The constants a„ and /?„ can be determined from the boundary conditions in Eq. (|28p . while the An are 
normalization coefficients. Using Eqs. ([^ . ([HO)) , and the orthonormahty of the ipn, Eqs. and (pS)) 
become: 



d oo 

-Pn(i) = - ^ i?nmA,„Pm(i) + q„(i) . 

?n— 



9t ' (33) 



where 



.L/2 .L/2 

Hnrn^ ds ds' "(An (s)Aiavg (s - s') V'm (s') ■ (34) 

J-L/2 J-L/2 

The matrix elements i?„„i can be evaluated through numerical integration. HWR neglect the off-diagonal 
portion of this interaction matrix since the diagonal elements typically dominate. However, as we will 
show later, this approximation leads to an inaccurate description of the simulation results for the end- 
monomer dynamics at short times, demonstrating that the off-diagonal elements are negligible only at times 
longer than the bending relaxation time. A more accurate approach is to take the whole matrix H, keep only 
the leading N x N sub-block (describing the interactions among the N slowest-relaxing modes), and exactly 
solve the resulting finite-dimensional version of Eq. (|33p . An appropriate value for N can be estimated as 
follows. For the oscillation described by mode V'n(s) from Eq. ([3T|) . the distance between successive nodes 
is approximately L/n. The high-frequency cutoff of this distance is on the order of two monomer diameters 
4a, so that only modes with n < L/Aa should be considered. Thus the natural choice is = L/Aa. In 
the results described in Sec. [5l we use this choice for all chains with L < 1600a. For longer chains with 
L > 1600a, calculation of the full matrix becomes numerically unfeasible due to roundoff errors in the highly 
oscillatory integrals of Eq. Thus for these chains we truncate A^ at the maximum value of N = 400. 

This approach gives accurate results at time scales much larger than the relaxation time of the n — 400 
mode, which is always the case for the time ranges of interest. 

To implement this approach, let J be the N x N matrix with elements Jnm = Hnm^m, A„ the eigenvalues 
of J, and G the matrix diagonalizing J: (CJC~^)„m = A„5„m. Assuming the eigenvalues A„ are distinct, 
and using the fact that H is symmetric, it can also be shown that the matrix C diagonalizes H through the 
congruent transformation: {CHC'^)nm = ©n^nm, defining diagonal elements Qn-^ If we introduce a new set 
of orthogonal functions ^'„(s) and the associated amplitudes Pn{t), Qn{t), 

N-l 



rn—0 
N-l 

Pn(0 ~ ^ ^ CnmPmit) , Qn(^) = ^ ^ Cnm^m{t) j 



(35) 



m— m— 



then Eq. ([33|) becomes 



^P„(t) = -A„P„(i) + Q„(t), ^^^^ 



{QrLi{t)Qmj{t')) = 2kBTSijS{t - t')&„Snm ■ 

This equation can be solved directly to yield the end-monomer MSD: 

(r2(t))^((r(L/2,<)-r(L/2,0))^) 

= 6i?<+^ A„(l-e-*/-"), 

n=l 

where the diffusion constant D = kBT<do'^o{L/2), the relaxation times t„ = A~^, and A„ = dkuTrnQn'^ni^/'^)- 
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Figure 3: The top panels show simulation and MFT results for the end- monomer and center-of-mass MSD 
{r'^{t))/a'^ for chains of length L — 100a and persistence lengths Ip = 10a, 20a, and 40a (a = monomer 
radius). The bottom panels show the local slope a = d\og{r'^ {t)) / dlogt for the end-monomer MSD, with a 
horizontal line at 2/3 marking the Zimm theory prediction. Two times are indicated by dotted vertical lines: 
Tp, the time at which the end- monomer {t'^{tp)) = Ip] and ri, the longest relaxation time of the polymer, as 
determined from the numerical simulations. 

5 Results 

The Brownian dynamics simulation and MFT results for the end-monomer and center-of-mass MSD are 
shown in Figs. [3][5]for chain lengths oi L — 100a, 200a, and 400a respectively, at various persistence lengths 
Ip. We also show in the bottom panels of each figure the effective local exponent a{t) = d\og{r^ (t)) / dlogt 
of the end-monomer MSD curve. 

We find in both the simulation and MFT results that a{t) passes through a minimum in the intermediate 
time range where Ip < {r'^{t)) < r1^{L). The location of this minimum is on the order of ti, the longest 
relaxation time of the polymer. For i ^ ri, as the end-monomer curve approaches the center-of-mass MSD, 
(r^jjj {i)) — 6Dt, the local slope a{t) tends toward the limiting value of 1. On the other hand, for t < Tp, 
where {r'^{t)) < Ip, the stiffness of the polymer dominates, and a{t) varies in the range w 0.8 — 0.9. We will 
discuss both the intermediate and the short-time regimes in more detail below. 

There is very good agreement between MFT end-monomer MSD predictions and simulation results in 
time regimes where simulation results have sufficiently converged to make a comparison. (For {r'^{t)) this 
comparison is possible for nearly the whole simulation time range; for a{t) the numerical uncertainty at the 
largest times becomes significant, and is on the order of the scatter in the plotted data points.) Additionally, 
dynamical parameters like the diffusion constant D and relaxation time ti determined from the simulation 
data compare favorably with their MFT values, as shown in Fig. [S] The values of D were obtained from the 
simulations by fitting the center-of-mass MSD data to the straight-line form (r^,j-j (t)) — 6Dt. To extract 
Ti, the autocorrelation function of the end-to-end vector was calculated, Ccc{t) — (rco(i) ■ rco(O)), where 
roo(i) = rA/(i) — ri(i)- For sufficiently large t, this function takes the form of a simple exponential decay, 
Cee(t) ~ exp(— t/ri), from which n can be estimated. For both D and n, only data points for which 
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Ip — 15a Ip — 30a 
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tkBTfj,Q/a^ 

Figure 4: Same as in Fig. [3l but for chains of length L ~ 200a and persistence lengths Ip ~ 15a and 30a. 
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Figure 5: Same as in Fig. [31 but for chains of length L = 400a and persistence lengths Ip = 30a and 60a. 
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Figure 6: Diffusion constant D = D / [ksT ^t)) (top) and longest relaxation time fi = tiUbT hq/ a? (bottom) 
as a function of Ip/ L for three different chain lengths L. The squares are Brownian dynamics simulation 
results, while the solid curves are calculated from MFT. 

convergence was achieved were included in the fitting (the criterion for convergence was that the local slope 
dlog(rL„.(i)>/rflogi«l.) 

The main discrepancies between the two approaches are in the local slopes of the MSD curves at the 
shortest times, t < 10^ a'^ /{kBTno). This can be explained by the fact that the small-scale motions at 
short times are particularly sensitive to the discrete nature of the polymer chain and the more strongly 
fixed monomer-monomer separation in the simulation, thus giving rise to differences with the continuum 
mean-field approximation. In fact we can make the MFT mimic the simulation more closely if we exclude 
the contributions of a fraction of the highest modes in the sum of Eq. (|37|) . by changing the upper limit 
from iV — 1 to A^' — 1, where N' — cL/Aa, < c < 1. A value of c w 1/2 gives the closest approximation to 
the simulation MSD and a{t) curves, irrespective of L and Ip. This roughly corresponds to excluding modes 
where the distance between nodes is shorter than four monomer diameters. The results are shown in Fig. [7] 
for three different chains, with the modified MFT labeled as MFT'. The long-time behavior is unaffected 
by removing the highest modes, but at short times the MFT' a{t) curves fit the simulation data much more 
closely. We rationalize this as being due to an effective cutoff of fluctuations at small wavelengths due to 
the spring stiffness in the simulation, which is not represented well by the Gaussian MFT elastic energy in 
Eq. (ini). 

In Fig. [7] we also show the MSD and a{t) curves calculated using the HWR model. This model follows 
the same basic approach as in Sec. 01 with two additional approximations: (i) only the diagonal elements 
Hnn of Eq. (l34l) are used; (ii) for n > 1, the i/„„ are evaluated approximately as^ 

„ fefioa L-s ( 3dM 

tinn~^\\ T] as^==exp cos a„s. (38) 

^ L vMs) V 2cr(s)y 

The net effect of these approximations is negligible only for i > ti, where the HWR and MFT results 
overlap. For t < ti there are significant differences with respect to the simulations. Here the HWR model 
overestimates the end- monomer MSD and underestimates a{t). The discrepancy is only slightly reduced 
by avoiding the approximation of Eq. (|38p : the main weakness of the HWR model is that the off-diagonal 
matrix elements Hnm are not included in the calculation. Taking these into account, as the MFT results 
demonstrate, gives a much more accurate description of the simulation data at short and intermediate times. 
Despite these differences, the sub-Zimm scaling regime exists in the HWR model, and the deviation below 
2/3 is even larger than in the MFT. (The sub-Zimm scaling is also implicitly evident in related quantities 
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Figure 7: Comparison of the simulation results for three semiflexible chains with several analytical ap- 
proaches: the heuristic scaling argument described in Sec. O the MFT described in Sec. HI and the MPT', 
HWR® models described in Sec.[5l The top panels show the end- monomer and center-of-mass MSD (r^(i)); 
the middle panels show the relative difference between the various theoretical results for the end-monomer 
MSD and the simulation data; the bottom panels show the local slope a = d\og{r'^ (t)) / dlogt for the end- 
monomer MSD. Two times are indicated by dotted vertical lines: Tp, the time at which the end-monomer 
(r^(Tp)) — Ip', and ti, the longest relaxation time of the polymer. The ti and Tp values shown are from the 
numerical simulations. 



14 




Figure 8: Results of the MFT for various semiflexible chains, (a) The end-monomer MSD effective exponent 
a = dlog{r^ {t))/d log t versus t/ri for polymers of diameter 2a, contour length L — 12000a, and various 
persistence lengths Ip = 50a — 3000a. The dot along each a{t) curve marks a{Tp). The time Tp is where 
{r'^{Tp)) = Ip, and ti is the longest relaxation time of the polymer, (b) The minimum value of a over all t 
versus Ip/L for several chain lengths L = 100a — 12000a. (c) The minimum value of a over all t versus L/lp 
for several Ip = 50a — 1800a. (d) The circles show the minimum value of a in each of the curves in panel 
(b), plotted as a function of 1/ log(i/a). Superimposed is a straight-line fit to the data points. The squares, 
with the dashed straight- line fit, are the corresponding results of the heuristic scaling argument, taken from 
Fig.Hd). 
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Figure 9: Comparison of the MFT and heuristic scaling results to experimental MSD data taken from 
Ref. 2 for three lengths of double-stranded DNA: (a) 2400 bp; (b) 6700 bp; (c) 23100 bp. The dimensional 
parameters in the theories are fixed at: a = 1 nm, rise per bp = 0.34 nm, Ip = 50 nm, T = 293K, viscosity 
of water at 293K ~ 1 mPa s. The top panels show the end- monomer MSD (r^(t)), and the bottom panels 
show the corresponding effective local exponent a = (ilog(r^(t))/dlogt. 
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calculated from the HWR approach, like the fluorescence correlation function studied in Rcf. 3.) 

For comparison, Fig. [7] also shows the (r'^it)) and a{t) curves calculated from the heuristic scaling 
argument of Sec. [H Despite its simplicity, it is able to capture the trends of the simulation and MFT data 
quite well, though for shorter chain lengths it gives a shallower dip in a within the intermediate regime. 

Given the success of the MFT at reproducing the simulation results, it is interesting to see what the 
theory predicts for longer chain lengths where Brownian hydrodynamics simulations become impractical. 
Fig. [Ha) shows a{t) curves for L = 12000a, Ip = 50a — 3000a, with the point a{Tp) on each curve marked by 
a dot. As in the shorter chains, there is a broad dip in a(t) between Tp and ti, but the minimum of a{t) has 
been shifted to below 2/3. In fact the dependence of this minimum on L and Ip, illustrated in Fig. [5fb-d), is 
qualitatively the same as that derived from the heuristic scaling argument in Fig.[IJb-d): there is a general 
trend of mint decreasing with L, and in particular the smallest value possible at a given L, min; mint ct, 
has a nearly linear dependence on 1/ log(L/a) (the heuristic result from Fig.[ljd) is also drawn for reference). 

At Ip — 50a, corresponding to the persistence length of DNA, mint ct ranges from 0.698 aX L — 100a 
to 0.617 at L = 12000a. Within this range we can make a detailed comparison for three particular chain 
lengths where experimental MSD data for double-stranded DNA is available from Ref. 2: L = 816a, 2278a, 
and 7854a, or equivalently 2400 bp, 6700 bp, and 23100 bp (using a = 1 nni, and a rise per base pair of 
0.34 nm). The experimental end- monomer MSD for these three cases is shown in the top panels of Fig. [HI 
The bottom panels show the local slope a{t), which can be estimated at each t by fitting straight lines to 
the log-log plot of MSD data points with times ti within a small range around t, defined by the condition 
\\ogiQ ti/t\ < 0.15. Together with the experimental results for the MSD and a(t) are the curves predicted 
by the MFT and heuristic scaling argument. Besides the length scale parameters mentioned above, the 
other dimensional variables in the system are set at the following values (taken from the literature and 
the experimental conditions): Ip — 50 nm, T — 293K, viscosity of water at 293K = 1 mPa s. At longer 
times (> 10 ms) there is quite good agreement between both theories and experimental data, particularly 
notable since there is no fitting parameter involved in the MFT. The discrepancies arise in the short and 
intermediate time regime, where the experimental MSD is consistently higher than the theoretical one, the 
difference increasing to roughly a factor of 2 — 3 at the shortest times measured. The discrepancy in the MSD 
is on the order of 0.001 — 0.01 /im^, corresponding to length scales roughly 30 — 100 nm. Experimental data 
seems to indicate faster displacement of the end monomer at very short times followed by slower increase of 
the MSD at intermediate times compared to the theoretical predictions. 

The effect of the higher experimental MSD is to push the local slope down relative to the theoretical 
value. Thus the intermediate dynamical regime, in the range 0.01 — 10 ms, is characterized by a broad region 
with a{t) close to 0.5, in contrast to the MFT results where mint en is between 0.633 for L — 2400 bp and 
0.617 for L = 23100 bp. Though there are large uncertainties in the experimental data for t < 0.01 ms (on 
average 50% for L = 2400 bp, going down to 10% for L = 23100 bp), the rough trend in the local slope 
appears to show a rapid increase in a{t) as t is decreased. This rapid crossover again contrasts with the MFT 
curve, where the increase in the local slope is more gradual. The heuristic scaling results support the MFT: 
with the crossover exponents and fitting constants set at the values shown below Eq. the heuristic {r'^{t)) 
almost perfectly overlaps with the MFT curve in all three cases, and the local slopes are consequently also 
very similar. Independently, we also checked if it was possible to find an alternative set of fitting parameters 
which would make the heuristic (r^(t)) agree with the experimental data, but we were unable to obtain a 
reasonable fit. 

Although the existence of an intermediate regime with sub-Zimm scaling is found in both experiment 
and theory, the quantitative discrepancy of the scaling behavior points to a gap between the experimental 
system and the theoretical approaches. Possibly, a semiflcxible polymer model based on a worm-like chain 
is suflacient only for describing the large-scale motions of the DNA. There may be some missing elements 
in the theory (for example an additional degree of freedom present in DNA, like torsional dynamics) that 
lead to faster motion at shorter scales. On the other hand, there is also the possibility that limitations 
in the setup and analysis of FCS measurements could contribute to the discrepancy. Deviations from the 
assumed Gaussian profile of the confocal detection volume and uncertainties in the diffusion coefficient of 
the rhodamine molecule used to calibrate the shape of this volume have been shown by alternative methods 
like two-focus FCS to lead to substantial systematic errors in the single-focus setup. The uncertainties in 
the FCS analysis are highlighted by the differing results produced by independent studies of similar double- 
stranded DNA systems: one yielding a substantial sub-Zimm regime,^ with local exponents near the Rouse 
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Figure 10: The values of (a) fvi = rnfesT/io/a^, (b) A„ = A„/a^ as a function of mode number n for 
a semiflexible polymer of length L = 12000a with various Ip = 50a — 3000a, calculated using MFT. The 
parameters A„ and relaxation times t„ determine the subdiffusive behavior of the end-monomer MSD 
through Eq. ([37|) . The insets show the effective local exponents describing the scaling of r„ and A„ with n: 
7n = -dlogTn/dlogn, 6n = -dlog An/dlogn. 



limit, and others giving a smaller deviation below the Zimm value over shorter time ranges,^' ^ as fitted by the 
HWR model. Regardless of these issues, there is one aspect in which both the experimental and theoretical 
approaches agree: an intermediate dynamical regime is present in the end-monomer MSD results, and this 
regime shows sub-Zimm scaling for long enough chains. 



6 Discussion 

To understand the intermediate dynamical regime in more detail, and to see where the deviations from the 
Zimm model arise, let us analyze the behavior of Eq. ([57)1 . whose general form is shared by any theory that 
expresses the end-monomer MSD in terms of contributions from normal modes. Both r„ and A„ decrease 
approximately as power laws in n, and we can describe this decrease through effective local exponents 7„, 
5n > 1: 

d\ogT„ dlogA„ 

In = --7j , Sn = Jj . 39) 

alogn dlogn 

Several examples of r„ and A„ are plotted in Fig. [10] for L = 12000a and Ip = 50a — 3000a. The variation 
of 7„ and Sn with n is shown in the insets. The behavior of the MSD slope a{t) can be directly related to 
these local exponents. Let us assume that within some time range r„2 <C i ^ t„j, n2 S> jt-i, the associated 
normal mode exponents are approximately constant: 7„ « 7, 6n ~ 6 for ni < n < n2. Then the dominant 
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contribution to the subdifFusive behavior at times t within this range is given by: 



n—ni 

^A^^nl dnn-*(l-e"*""/""i"^) (40) 



A, 



In the second hne we extended the upper hmit of the integration from n2 to oo using the fact that t ^ t„2 
n2 ^ ni, and in the third hne E^{x) denotes the exponential integral function, E^{z) = dt exp{— zt) / 1" 
For t <C T„^ , Eq. can be expanded to the leading order as 



2u\\ _ ^ni"-lT. /^l -(5 + 7 



This imphes that for t„2 ^ i <C the local slope is given by a{t) « ((5 — l)/7. As can be seen in the insets 
of Fig. [ini there are two distinct regimes for 7„ and (5„: one for modes n ^ L/lp, corresponding to length 
scales greater than Ip, and another for modes n ^ L/lp, corresponding to length scales smaller than Ip. A 
continuous crossover occurs from one regime to the other for n ~ 0{L/lp). These two regimes in turn lead 
to differing behaviors for a(t). We will consider each regime separately, focusing on earlier predictions for 
each case and how they compare to the present results. 

Modes with n <C L/lp correspond to internal polymer dynamics on length scales between Ip and L, and 
this is precisely the intermediate dynamical regime that we have mentioned earlier. In the simplest analysis, 
ignoring hydrodynamical effects, these modes should be described by the Rouse model, particularly in the 
flexible limit of small n where the length scales are much greater than Ip. The Rouse theory yields the 
following expressions for t„ and A„,^ 

^»°"^°- 1o 2, t 2 "'^ Af- = ^il^n-^ (42) 



where b is the Kuhn length, b — y^r'^^{L)/M. The local exponents in the Rouse model are constants: (5„ = 2 
and 7„ = 2. Using Eq. (j4T|) . we find the following asymptotic behavior for the end-monomer MSD: 



{r^{t)) » Af™=''T(l/2)(t/Tf°""<=)i/2 

t 



m'kBT^,oV^\,/, (43) 



TT 



This is the origin of the Rouse scaling result a(t) — 1/2. In the presence of hydrodynamic interactions, the 
Zimm model is expected to hold, with Eq. (|42p modified as:^ 

^Zimm ^ b^{L/7:a)'^^^ n^^/'^ ^Zimm _ ^Rousc ^^^-j 

" 12y/6kBTiiQa " " 

Here the local exponent 7„ is 3/2, and thus the asymptotic behavior of the MSD becomes: 

^^2^^^^ ^ 12r(l/3)(2afcgT/.o)^/3 ^,^3 ^ ^^^^ 

leading to the Zimm scaling a{t) ~ 2/3. 

The MFT calculations, however, give a different picture, deviating from the Zimm result. Consider 
the chain in Fig. [TO] closest to the fiexible limit: L = 12000a, Ip = 50a. The exponents 7„ and (5„ are 
approximately constant for n < 10, but are shifted from the Zimm values: 7 ~ 1.74, 6 ~ 2.04 averaged over 
n — 1 — 10, giving a = {d—l)/j ~ 0.60. Indeed in the corresponding local slope curve plotted in Fig.lHJa) the 
a{t) value is nearly constant over the time scales associated with these modes {t/n « 0.018 — 1), reaching 
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a minimum of 0.617. It is these shifts in 7„ and Sn from the Zimm theory predictions that lead to an 
intermediate dynamical regime for longer chains where a(t) < 2/3. 

To get an analytical estimate for these shifts within the framework of the MFT theory, one can approxi- 
mately evaluate the integrals for the interaction matrix elements Hnm in Eq. (|34p , and account for the effects 
of the off-diagonal elements using perturbation theory. The details of the approximation can be found in 
Appendix B. For n <C L/lp the results are: 

3pi{n) + 5.07p2{n)J^K (^) 2pr,{n) + 0. 0393^6 Ha/J A' (^) 

^ ^ " ^ V= (46) 



2p3{n) + 2.89p4{n)^K (^) pr{n) + 0.0l57ps{n)^K 



where K{x) = y ^ (^1(2;) — ^i(3/2)) — 4y + 2, and pi{n), . . . ,psin) are polynomials in n of the form 

Pi{n) = 1 + tti/n + bi/n^ with coefficients and bi given in Table 1 of Appendix B. For L ^ a and 

Ip approaching the flexible limit, Ip — > 2a, the first terms in the numerators and denominators of the 7„ and 
Sn expressions dominate, and thus there is a range of modes 1 <C n <C L/lp where 7„ « 3/2 and 6n ~ 2, 
in agreement with the expected Zimm scaling for a flexible chain. However, for a semiflexible chain where 
Ip ^ a, corrections to the Zimm values become more important. Using the fact that Ei{z) w — logz as 
z — > 0, the K{6a'^/lp) terms in the 7„ expression lead to a positive shift of order (nlp/Ly/"^ \og{lp/a). For Sn 
the shift upward is smaller, of order [Ip/nLY/'^ \og{lp/a). These corrections due to semiflexibility are evident 
in the exact numerical results for chains of length L = 12000a shown in the insets of Fig. [101 particularly for 
Ip = 50a and 100a where a n ^ L/lp regime is identifiable. As expected from the analytical approximation, 
the deviation in 7„ from the Zimm value is more significant than that of Sn- In fact the averages of 7„ and 
Sn for n = 1 — 10 from the approximate expressions in Eq. (|46p are 1.68 and 2.06 respectively, comparable 
to the numerical results 1.74 and 2.04 quoted above. 

In the other regime, for modes with n ^ L/lp, the oscillations are at length scales smaller than Ip, where 
the rigidity of the chain is the dominating factor. For this case it is easiest to consider first the MFT in the 
absence of hydrodynamic interactions, and then see how the final results are modified when the interactions 
are included. In the free-draining limit, the interaction matrix Hnm = 2a/xo5„m, and for large n the constants 
a„ in Eqs. ([^ - ([5^ are approximately a„ ~ n-n/L? With these simplifications we find 

n-\ A„«i^n-\ (47) 



SajiolpkBTir'^ ' Ipir^ 
for n » L/lp. Thus Sn — ■^n — 4. Plugging these results into Eq. gives 

for the end-monomer MSD. The scaling {r^{t)) ~ <3/4 

is a well-known property of monomer motion in the 
stiff-rod limit, as seen in theory, simulations,^^' '^"■'^^ and experiments.^' '^^"'^^ Though hydrodynamic 
effects are typically expected to induce only weak logarithmic corrections in this limit, we find that including 
these effects in the MFT does have an observable consequence. In the insets of Fig. \W\ all the 7„ and Sn 
curves appear to overlap for n 3> L/lp, but their values are shifted away from 4: 7„ gradually decreases with 
n, varying between 3.9 and 3.5 in the range shown, and 5„ « 4.15 — 4.2. The behavior of 7„ and Sn lead to 
a{t) > 3/4 in this regime, as is seen most clearly in the large Ip results in Fig. [H^a), which exhibit a broad 
region where the a{t) curves converge over the range 0.8 — 0.85 for t <C Tp. 

The MFT analytical estimate for the n 3> L/lp case, using the approximation detailed in Appendix B, 
gives: 

12 24 

^""'^^-121og(^)' " ' + (3 + A - 12 log {^)){5 + A- 12 log (^)) ' ^^'^ 

where the constant A = I87 - 2\/67r + 6E1 (|) + 61og(6) « 13.06 and 7 w 0.5772 is Euler's constant. 
The functional forms for 7„ and Sn in Eq. (j49p . independent of Ip, describe the curves toward which all 
the 7„ and Sn results in the insets of Fig. [TO] converge for sufficiently large n, with 7„ shifted below 4 and 
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Sn shifted above 4. The gradual decrease in 7„ with n is similar to an earlier theoretical approach where 
hydrodynamics was explicitly considered: in Ref. 11 the relaxation times for a stiff-rod were found to scale 
like T„ cx n~^/log(L/an7r), corresponding to 7„ = 4 — 1/ log(i/an7r). 

In a more intuitive fashion, the heuristic scaling developed in Section II allows to trace back the deviations 
from the traditional Zimm and worm-like-chain scaling results in Eqs. (HS)) and to the slow crossovers 
in the diffusion and spatial size of sub-chain segments. 

7 Conclusion 

Between the flexible and stiff-rod limits hydrodynamic interactions modify the scaling of the end-monomer 
MSD in ways that are not accounted for in the Zimm model, or in earlier semiflexible polymer theories. In 
particular, there exists an intermediate dynamical regime for sufficiently long polymers with local exponent 
a{t) between 2/3 and 1/2, the Zimm and Rouse predictions. We have investigated this regime through a 
worm-like chain model, in conjunction with a variety of theoretical techniques: Brownian hydrodynamics 
simulations for shorter chain lengths, supplemented by mean-field theory with hydrodynamic pre-averaging 
for longer chains where the simulations are not practical. In the cases where both MFT and numerical 
results are available, there is very good quantitative agreement between them. The two approaches are 
further supported by a heuristic scaling argument that can accurately capture the trends in {r'^{t)) and a{t) 
and that allows us to connect the observed sub- Zimm scaling regime to previous scaling approaches developed 
for the stiff- rod and the flexible-chain limits. Note that previous less accurate mean-field approaches that 
were used to analyze the FCS data of Ref. 5 give a sub-Zimm scaling range even more pronounced than 
found by us. 

Though the MFT and heuristics show a noticeable dip below the Zimm exponent of 2/3 at intermediate 
times, they do not reach the Rouse-like value of 1/2 seen in the experimental double-stranded DNA results 
of Ref. 2. Comparison between the experimental data and the theory raises a important issue: while the 
long-time data, corresponding to the large-scale dynamics of the DNA, is described surprisingly well by 
the MFT, serious discrepancies arise at shorter times. The small-scale motions revealed by experiment are 
significantly faster than predicted, indicating either a deficiency in the simple worm-like chain description 
or in the analysis of the FCS measurements. Further work is thus necessary in order to gain a complete 
understanding of the monomer dynamics of DNA in solution. 
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Appendix A: Mean-Field Theory of an Extensible Worm-Like Chain 

As an alternative to the mean-field theory of Sec.[4l which begins with the inextensible Kratky-Porod chain of 
Eq. (|12p , we can derive a mean-field model based on the extensible worm- like chain Hamiltonian used in the 
Brownian dynamics simulations, Eq. (jlOp - (llip . thus making explicit the relationship between the simulation 
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and analytical results. Ignoring the Lennard- Jones term, the simulation Hamiltonian has the form 



where d — 2a, T — ^d? , and u; = (r^+i — Yi)/d. For large F (the case in the simulations), the values of 
Ui = \u.i\ « 1, and we can expand [ui - if = (yr+~(uf~T) - 1)^ w (uf _ 1)2/4 _^ 0((m| - 1)^). Keeping 
the leading term, we rewrite Eq. (jSOp as 

p A/-1 A/f-l 

^ ^ M E - 1)' - ^ E • + ^(*^ - 2) • (51) 

i=l i=2 

Neglecting the last term of Eq. ((5T|) . since it is a constant, the chain partition function Z is given by 

A/-1 

n dn^^-P"" 

(52) 



We can rewrite the integrand of Z using the relations 



e-|5("?-i)^ oc r d\ e-'3^-^(»^i)+^^? , (53) 



— zoo 



and Wi ■ Ui_i = i (m^ + — (u^ — Uj_i)^), where we have introduced an auxiliary variable for each i. 
The result, up to a constant prcfactor, is 

nlOO M-l 

Z^ \{d\,e-^^^^^^}\ (54) 



with 



and 



» Af-l 

F{{\}) = log j n dn, e-^^({^-» , (55) 



A/-2 A/-1 

^2 



c/({A.}) =d - ^) + i E ("^ - 

'i=2 i=2 



2d' 



i=l i=l 

To derive a mean-field model we can now apply a stationary phase approximation analogous to the one used 



in Sec. 1123,24 



pioo M — 1 

Z^ II dX, e-'^^"^'" « e-^^«^"» , (57) 



where {Xf} satisfy 

dF 
dX, 



= 0, i = l,...,M-l. (58) 
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From symmetry, we know A^' must have the property A^' = ^Hi-i for ^^^d thus we can write the solution 
to Eq. ([nH) in the form 

^f = ''^ + l2' i^2,...,M~~2, 

u , (59) 

for some set of values {i^i} and z/q, where i>i — vu-i- This yields a mean-field free energy 

Fmf = F{{\f}) - log / J] du, e-/^^"- , (60) 

i=i 

with 

J\/-2 

J/mF =rf ^ i^jMi + X! ~ Ui-l)^ 



2d 

1=2 j=2 



A/-2 



+ j/o(u? + Um_i) - d ^ - 2z^o (61) 



i=2 



- — 2. i'^^ + rfjj - — It + m^J • 

1=2 

In the continuum limit d — > 0, M — s- cx), Afd — s- L and we replace u^, Ui by continuous functions u(s), v{s) 
of the contour variable s. Assuming the Hamiltonian parameters e and F remain fixed in this limit, we have 



Fmf = - /3~^ log j Vue-f^^^^ , 
Umv ^ I dsv{s)u^{s) + - i ds 



du{s) 

-L/2 2 y_i/2 V 

rL/2 (62) 

+ j/o(u2(i/2) + u2(-L/2)) - / dsiy{s) 

J-L/2 

-2vo - — I ds v{s) C , 

i J-L/2 i 

where C = limc(^o(2Le2/r(i2 + e^/Fd) is an infinite constant independent of z/(s) and vq. The stationary 
point condition Eq. ([58]) becomes 

^—-=0, — = 0. 63 

Physically Eq. (l63|) implies the following constraints: 

(u2(s)> = l + ^, -L/2<s<L/2, 

2. (64) 
(u2(-L/2)> = (u2(L/2)> = l + -, 

where ( • ) denotes the thermal average with respect to the Hamiltonian L^mf- As expected, the magnitude 
of the tangent vector fluctuations become smaller as the extensibility parameter F increases, going to the 
limit {v?{s)) — 1 for all s when F ^ oo. As will be seen below, in this limit the present theory reproduces 
the results of Sec. SI 

Unfortunately Eq. is not analytically tractable, because the Fmf given by Eq. ([S^ cannot be 

evaluated in closed form for an arbitrary function i'{s). On the other hand, the stationary point condition in 
the discrete system, Eq. ([58|) for F in Eq. ([55|l . can be solved numerically for small M, and a representative 
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0.12 



0.11 



0.10 




Figure 11: Numerical solution to the stationary point condition Eq. ([58]) for a system with M — 20, 
e = 10 ksTd, and T = 400 ksTd. The solution is expressed in terms of Ui = Xf — e/cP for z = 2, . . . , M — 2, 

_^-e/2d? is 0.7831 /csT/d. 



in units of ksT/d. The value of Vf)/d = \f - t/2d^ = Xf,j_ 



set of i^i are shown in Fig. [TT] (for M = 20, e = lOfc^Td, F = AOOkBTd). For large F > e, the are nearly 
constant for 2 < z < M — 2, and we can use this fact to make the following approximation in the continuum 
limit: replace i/(s) by a constant ly in the Hamiltonian Umf of Eq. ((5^ . Thus the first part of Eq. ([55]) 
becomes dFuF/di' = 0, implying a global constraint 



^^^d.(u^(s)) = (^l + -)L 



(65) 



instead of the local constraint (u^(s)) = 1 + 4e/F in Eq. ((64l) . 

This approximation, which becomes exact when F ^ oo, allows us to find a closed form expression for 
Fmf- Using a mapping of the first two terms of Umf to the quantum mechanical harmonic oscillator (with 
mass m — e and frequency uj = j e)^'^ the path integral for i^MF can be evaluated, giving the free energy 



MF 



-Lv — 2vf) 




log 




(3^/vecsc)\ I L\ — 



(66) 



/3-^j/oV2i'ecoth ( L\j — 



up to an additive constant. Using Eq. the stationary point condition dF^F/d^ — dFuF/di^o = can 

be solved numerically for v and vq given L, e, and F. For L ^ e, the condition takes the simple limiting 
form. 



— = —kb-i 
2 4 



F 



3 r(r + 6e) 

4 ^ (F + 2e)(F + 4e) 



(67) 



When F ^ oo, Eq. l[67|l reduces to Eq. (jlSp in Sec. [H and this is generally true of the stationary point 
condition for any L. 

The equilibrium properties of the chain described by the Hamiltonian Umf can be calculated with a 
similar approach to the one used in Ref. 27, where distribution functions were derived for the F = cx3 mean- 
field model. The main result is G(s, s'; x; u, u'), the probability density for finding two points on the chain 
at s and s' with s' > s, having spatial separation r(s') — r(s) ~ x, and tangent vectors u(s) = u, u(s') = u'. 
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The full expression for this probability is 

G{s, s'; x; u, u') = 

/ Ais,s') 
V47r3B4(s' - s) 



exp 



^^(^ 5^(737^ 

C(s + L/2)u2 - C(L/2 - s')u' 



where the functions A[s), Bi{s), i = 1, . . . ,4, and C(s) are given by 

A(s, s') = 4(Bi(s' - s) + C(s + L/2)) 



•(SiGs'-.s) + C(L/2-s'))-S|Gs'-s), 
i?i(s) = coth(sw) , -82(5) = /3ew csch(sa;) , 

s euj /suj\ 
(3eLo{euj + 2(^0 coth(sw)) 



41^0 + 2euj coth(sa;) 



with LO = yj2v/e. 

From 0(5, s'; x; u, u') we can calculate other properties of the chain, for example the tangent vector 
correlation function (u(s) • u(s')) for s' > s, 

(u(s) • u(s')) = f d^xu(s) • u(s') Gis, s'; x; u(s), u(s')) 
For L ^ e, —L/2 <C s, s' <C i/2, Eq. ([70)1 can be simplified using Eq. ([57)1 for 1/ and fo, giving 

When r ^ CX3, the tangent correlation function reduces to the Kratky-Porod form of Eq. (|19p . (u(s) •u(s')) = 
exp(-3(s' - s)kBT/2e) = exp(-(s' - s)/lp). 

Thus we have shown that a mean-field theory based on the extensible worm-like chain Hamiltonian used 
in the simulations gives results very similar to the MFT described in Sec. IH with the finite extensibility 
leading to small corrections to the parameters of Umf on the order of e/F. 

Appendix B: Analytical Approximation for 7„ and 5n 

In order to derive analytical expressions for the exponents 7„ and (5„ from the MFT theory, we make several 
approximations in the derivation described in Sec. 3) Since the off-diagonal elements of the interaction matrix 
Hnm defined by Eq. (I34p are smaller than the diagonal ones, we will treat them as a perturbation. To first 
order in the perturbation expansion, we can write the following expressions for r„, 8^, and ^'„(i/2) when 
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n>0: 



tt2 \ 

_ A-1 ^ \-lf7-l _ H^-2\-l \^ -"nm^ni 



9„^g„„ + 2^ Z^-^" , (72) 

_rz 7^7^;, An J I mm ■^m 

*„(L/2) = ^ V„(X/2)(C-1)™„ « Vn(V2) + + E ^rn(L/2) . 

V ^ ^7171 , ^7\ri^7\ ^7rim^rn 

m rn^n 

From these expressions one can also calculate — dkBTrnQn^ni^/'^)- The double integral for i^nm in 
Eq. (j34p can be rewritten in terms of new variables h — s — s' and — s + s' as follows: 

.TO. Z"^^" ^ Z"^^' ^ , w ^0(|s-«'|-2a) / 6a2 \ ^ ,^ 

= 2afiodn„i + \ -afiQ / dh / — - — ^— — exp — — tt 

Since <t(/i) = 2lph — 21^(1 — e~^/^p) can be approximated as (t(/i) « ft.^ for h <^lp and cr(ft.) « 2Zj,ft, for h ^ Ip, 
we can split up the h integral above into two pieces: 

Hnrn ~ 2a^o'5„m + a^o(4m + -^im)> (74) 



where 



2a \ J J-L+h \ 2 ) \ 2 

r(2) /6 /"^^^ 1 / 3a2\ /•^-'^ ^ , + , fw~h 

^n,n = \l - I dh —^== exp ( -— ) / dwipn y — ^ — j p. 



n Ji^ y^2lph V ^P>^ 



(75) 



To complete the approximation, we will estimate these integrals in the two mode regimes discussed in Sec.[6l 
one for the case n <C L/lp, the other for L/lp. Since the biggest perturbation contributions in Eq. ((72|) 
for r„, On, and '^n{L/2) come from states with ni in the vicinity of n, it is sufficient to consider matrix 
elements Hmn for ni in the same mode regime as n. 

B.l. n,m ^ L/lp regime 

In the limit of long chain lengths, where L 3> Zp, a, the functions ipnis) and constants a„, A„ for < rt <C L/lp 
in Eqs. ([ST]) and ([32]) simplify to: 

r(_l)("-i)/2 /|sin(i:g£) nodd, 

|^(_l)n/2^2^os(2f^) neven, ^^g^ 

Due to the symmetry of the ipnis) functions, the matrix elements Hnm are non-zero only when n and m 
are both odd or both even, so the perturbation expansions in Eq. ([72)1 can be done independently for even 
and odd states. For simplicity, we will assume n and m are odd for the rest of the derivation. Carrying out 
the analogous approximation for even n, m, will lead to qualitatively similar final expressions for and 7„, 
with slight shifts in the numerical coefficients. 
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Plugging Eq. (|75|) into Eq. (|75p . we can approximately evaluate the integrals in the large L limit: 

2ae-3/2 - Ipc-^^^'^l + a\/67f erf (1/372) - 0^67^ erf(V6a//p) 
'f [i?! {Qa'/lD-E, (3/2)] 




n ^ m, 
n = m. 



/(2 

nm 



2VI 



(77) 



Plugging these results into Eq. (|7^ for i/,„„, we can also estimate the sums involved in the perturbation 
expansion of Eq. (|72p : 



E 



tt2 \ 



/-1087r + 277r2 -126 + 727r - 16\/3- 



7T 



Hrin An Hmrn Am 



1 / /-\ -9 + (9 + \/3) TT 72 - 457r + 4V37r 

— (^-9 + 4V3J + ^^-5^^ 216^3^2 



(78) 



Combining the results of Eqs. ([72|) . ([74)) . and ([76)1 - ([78)1 . we can derive the following expressions for 7„ and 

5n- 



In 



6n 



d\og n 
dlogn 



n dTn 

T„ dn 



3pi(n) + 5.07p2(^) 



-K 



2p3{n) + 2.89p4{n)^K (^) 



dA„ 2p5(n)+0.0393p6(n)V;S:^n 



6a' 



(79) 



where 



^ (i?i(a;)-i?i(3/2))-4W- + 2, 

TT V TT 



(80) 



and pi{n), . . . ,psin) are polynomials in n of the form Pi{n) = 1 + ai/n + bi/n'^ + ■ ■ ■ . The first two coefficients 
tti and bi are given in the following table: 



i 


1 


2 


3 


4 


5 


6 


7 


8 


0.1 


-0.176 


-0.0919 


-0.297 


-0.148 


0.282 


0.608 


0.188 


0.435 




-0.00186 


-0.00504 


0.0179 


-0.00196 


0.0440 


0.185 


0.0220 


0.103 



Table 1: Coefficients a.^ and bi of polynomials Pi{n) — 1 + ai/n + hijv? -t- 



B.2. n,m^ L/lp regime 

With the assumptions that L ^ Ip and Ip 3> a, the functions tpnis) and constants a„, A„ for n ^ i//p in 
Eqs. dSI]) and ([^ become: 



Ipnis) 




\ 2L ; 



1 sinh(-n-(2rt-l)s/2L) 
L cosh(7r(2n-l)/4) 



n 



odd, 



/ 7r(2n-l)s\ , /l cosh(7r(2n-l)s/2L) 

i 2L j + V ^ sinh(^(2»-i)/4) "even, (g^) 

2,kBTlpTT^{2n~ If 
32I4 ■ 

Again we will focus for simplicity on the case of odd n and m. For n,m^ L/lp, the first integral in Eq. (|75p 
dominates, Inm ^ ^nm, so we can write Hnm ~ 2afj,oSnm + ap-olnm- The integral /^m can be approximated 
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as: 

2V6 (™''W)(log(iS2:^£5ii^)+27)-(m+n-l)(m+n)^7r 
j(l) ^ ^ (m+n- 1) (m+n) (m^ _|_„2 ) n ^ m, 



V6 



nn ■ 



{37 + ^1 (f) +21og(^%™)} 
isulting sums in the perturba 



n = m, 



where 7 « 0.5772 is Euler's constant. The resuhing sums in the perturbation expansion of Eq. (j72p are: 



^ HnnK-H^^X^ ~ ny^(3 + A- 121og(^)) 



. / , X (83) 



^ "^"^ / ' Hnn\n - Hmra^m ^ VI I 3 + A - 12 log (-^) i ' 

where the constant A — I87 — 2\/67r + 6i?i (|) + 61og(6) ~ 13.06. This yields the foUowing expressions for 
7„ and (5„: 

12 24 

^""'^^-121og(^)' ^""'^(3 + A-121og(^))(5 + A-121og(^))- ^''^ 
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